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Positronium spectroscopy (n = 1 hyperfine splitting, n = 2 fine structure, and the 1S-2S interval) 
has reached a precision of order 1 MHz. Vigorous ongoing efforts to improve the experimental results 
motivate the calculation of the positronium energy levels at order ma 7 . In this work we present the 
result for a complete class of such contributions-those involving virtual annihilation of positronium 
to three photons in an intermediate state. We find an energy shift of 2.6216(ll)ma 7 /(rnv) 3 for the 
spin-triplet S state with principal quantum number n. 

PACS numbers: 36.10.Dr, 12.20.Ds 


Positronium, the bound state of an electron and its 
antimatter partner, the positron, is of fundamental in¬ 
terest for a number of reasons. It represents the purest 
example of binding in quantum field theory as the con¬ 
stituents are structureless point particles and their low 
mass implies that the dynamics is dominated by QED- 
strong and weak interactions effects are negligible. As 
a particle-antiparticle bound state, positronium exists in 
eigenstates of the discrete symmetries parity and charge 
conjugation, leading to the possibility of leptonic tests of 
the corresponding symmetries and also to the presence 
of real and virtual annihilation into photons. An addi¬ 
tional difference between positronium and other bound 
states such as hydrogen and muonium is the fact that 
the m/M recoil effects in positronium are maximal due 
to the equality of electron and positron masses. Crucially, 
positronium is accessible to high-precision measurements 
of its spectrum, decay rates, and branching ratios, and 
so provides an ideal laboratory for testing fundamental 
theory and calculational methods for bound state QED. 
Confirmed discrepancies between theory and measure¬ 
ment could be a signal for new physics in the leptonic 
sector. 

Several measurements of the positronium n = 1 and 
n = 2 transition frequencies have been done with un¬ 
certainties of order 1 MHz. These include the n = 1 
hyperfine splitting (l 3 Si — TSo), the 1S-2S interval 
(2 3 S'i — l 3 <Si), and a number of purely n = 2 intervals 
[lj. For example, the n = 1 hyperfine (hfs) results with 
the highest precision are CHS] 

AKhf s (Brandeis) = 203 387.5(1.6) MHz 

AKhfe(Yale) = 203 389.10(0.74) MHz (1) 

AE hfs (Tokyo) = 203 394.2(1.6) stat (1.3) sys MHz. 

The best result for the 1S-2S interval is jBj 

A£(2 3 5i - l 3 *Si) = 1233 607 216.4(3.2) MHz. (2) 

There is currently a significant push to explore new and 
improved approaches to the hfs measurement IMS] and 


to improve the 1S-2S result [141117] . At the moment, the 
hfs measurement is the more precise, but there appears 
to be more potential for improvement in the 1S-2S mea¬ 
surement. The ~ 1 MHz uncertainty in the hfs number 
represents a line-splitting factor of better than one part 
in a thousand; if a similar splitting of the 1S-2S transi¬ 
tion could be achieved (natural line width ~ 1.3 MHz), 
an ultimate experimental precision of a few kHz might 
be in reach. Significant improvements in the n = 2 fine 
structure [I5H22] would also be possible if a similar line 
splitting could be performed. 

On the theoretical side, all energy contributions 
through terms of 0(ma e ) are known analytically [23ff26] . 
In addition, the leading log terms of 0(ma [ In 2 (1/a)) are 
known fTH - Elj] , as well as the subleading log contribution 
to the hfs [50H52] . and a number of pure 0(ma 7 ) con¬ 
tributions [55H38] , The present theoretical result for the 
hfs is consistent with the latest of the hfs experiments 
(“Tokyo” in Eq. 0) but is ~ 3 <t above the earlier re¬ 
sults. The theoretical result for the 1S-2S interval is in 
marginal agreement (2.4<r) with the experimental value. 

The known contributions to the positronium n = 1 en¬ 
ergy levels are shown in Table |T] listed by order in a and 
ln(l/a). Also shown are the numerical values, in MHz , 
of ma 4 , ma 5 ln(l/a), etc. One can see that the actual 
energy contributions are somewhat smaller, by a factor 
of ~ 3, than the values ma. x \n v (\/a ) of the associated 
order. However, pure ma 7 contributions as large as sev¬ 
eral tenths of a MHz have been shown to arise when 
“ultrasoft” photons (energy and momentum ~ma 2 ) are 
involved [33, 34]. Knowledge of the full ma 7 result is 
essential for the unambiguous interpretation of present, 
and certainly future, measurements of the positronium 
transition energies, especially for the 1S-2S interval. 

Positronium energy contributions can be classified by 
whether or not they involve a complete annihilation of 
the electron-positron pair into some number of virtual 
photons as an intermediate state. Processes involving 
annihilation into one, two, three, and four photons all 
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contribute at Oirna 7 ). Charge conjugation symmetry 
implies that the £ = 0 spin-triplet states n 3 Si are only 
affected by one- or three-photon intermediate states. 
The one-photon-annihilation contribution at 0(ma 7 ) is 
known [M|. The purpose of this letter is to present the 
result for the 0{ma 7 ) energy level correction due to all 
processes involving three-photon intermediate states. 

The main challenge in dealing with the three-photon- 
annihilation (3yA) graphs is the separation of real and 
imaginary parts. Our approach is to deal with the graphs 
non-covariantly by first integrating over the energy com¬ 
ponents of the annihilation photon loop momenta. We 
then identify the terms that contribute imaginary ener¬ 
gies by Cutkosky analysis [35] and isolate the imaginary 
parts as coming from ln(—1 — ie) = —in. We illustrate 
this procedure in a calculation of the lowest order 3yA 
graph shown in Fig. [l] The diagram shown represents 
one of the six permutations of internal photons that con¬ 
tributes at this order. 

The energy shift through 0(ma 7 ) due to 3yA graphs 
depends on the spin of the positronium state and the 
value of the wave function at spatial contact c/o-other 
features of the wave function do not enter. All internal 
momenta are “hard”-of 0(m)-while momenta of order 
mo. and mo? do not contribute. Consequently, any con¬ 
venient bound state formalism can be used. We employ 
the formalism of Ref. m, in which the energy shift is an 
expectation value 


A E = i'S/SK'S’. (3) 


The two-body positronium states are, in this approxima¬ 
tion, 


^ -t <£o 



vF -> 00 



( 4 ) 


TABLE I: Known contributions to the average IS energy and 
n = 1 hfs at various orders in the perturbative series. The first 
column shows the orders that contribute, starting at ma 4 . 
(The leading contribution to the average energy — mot 2 /A is 
not shown.) The second column gives the numerical value of 
ma 4 , ma 5 ln(l/a), etc. (The energies are expressed in MHz 
using the 2014 CODATA recommended values for R^c and 
a.) 


order 

value 

-Eavg (IS) 

E hfs (l 3 Si - 1 4 S 0 ) 

4 

ma 

350377 

38 322.493 

204 386.630 

ma 5 ln(l/a) 

12 580 

3 003.302 

0 

ma 5 

2 557 

-1018.784 

-1005.497 

ma 6 ln(l/a) 

91.8 

2.869 

19.125 

ma 6 

18.7 

3.000 

-7.330 

ma 7 ln 2 (l/a) 

3.30 

-1.091 

-0.918 

ma 7 ln(l/a) 

0.67 

- 

-0.323 

ma 7 

0.14 

- 

- 


where x is the two-by-two two-particle triplet spin matrix 
X = 3 ■ e/ 1 /2 with £ the orthopositronium polarization. 
The relative momentum vanishes in this approximation, 
and 0o = \/m 3 a 3 /(8nn 3 ) is the wave function at spatial 
contact for a state of principal quantum number n and 
orbital angular momentum £ = 0. 

The explicit expression for the energy contribution 
shown in Fig. |T]is 


A E = (-l)*0g 


d 4 p d 4 g -i -i -i 
(27t) 4 (27t) 4 p 2 (P — p — q ) 2 q 2 


x tr 


0 0 


x tr 




(—iey^ 1 ) 
(-*e 7 P2 ) 


x (— ie^ 2 ) 




7 (—P/2 + q) — m 
i 

y(—P/2 + q) - m 


7(P/2 — p) — 771 


( ie'y ^ 3 ) 


0 X 
0 0 


( 5 ) 


where P = (2m, 0) is the positronium 4-momentum in 
the center-of-mass frame, p and q are the 4-momenta 
of two of the virtual photons, and the initial (—1) is a 
fermionic minus sign. We see that A E is proportional to 
ma e times a dimensionless quantity obtained by scaling 
a factor of m out of each momentum. The energy shift 
can be written as 


A E = 


ma 6 2 
7T 2 3 


dpdq p 2 q 2 



dp 0 dq 0 Tig 

2ni 2ni D\ a 


( 6 ) 


where p and q now stand for the magnitudes of the 
dimensionless 3-vectors p and q., u = p • q, T\ a is a 
product of two traces and the denominator factor is 
Dia = (p 2 + *ei)((2n — p - q) 2 + ie 2 )(q 2 + ie 3 )((p - n) 2 - 
1 + ie 4 ) 2 ((g — n) 2 — 1 + *e 5 ) 2 where n = (1,0). We have 
used distinct values for the various positive infinitesimals 
ei to facilitate evaluation of the energy integrals, which 
were performed as contour integrals. We extended the 
Po and q o contours by appending infinite semicircles in 
either the upper or lower half planes on which the inte¬ 
grand vanished. The closed contour integrals were then 
evaluated by use of the residue theorem. 
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FIG. 1: The lowest order 3yA contribution. 
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After performing the po and qo integrals we obtained a 
number of terms-those giving rise to the imaginary part 
all contained a denominator factor p+g + s — 2 — ie, where 
s = \p + q | and e here is a positive infinitesimal formed 
from ei, £ 2 , etc. We isolated the singularity by writing 

_ 1 _ = P + q-s- 2 

p + q + s — 2 — ie —2pq(u — u — ie) 

where u is the special value of u for which p + q + s = 2 : 

B = -l+ 2(1 -P)( 1 -‘,> = l- 2 (P + ‘>- 1 > , (8) 

pq pq 


The quantity u has physical values —1 < u < 1 when 
p < 1, qr < 1, and p + q > 1. This region, the or- 
thopositronium “decay triangle”, corresponds to physi¬ 
cal (positive energy) photons with total energy equal to 
that of positronium (2m in this approximation). The u 
integral was written as 


[' du 

u — u — ie 


= f 1 du f(p,q,u) - fjp.q^u) 

J- 1 u-u 

/■' 1 

+ f(p,q,u) / du ---(9) 

J_ 1 u — u — ie 


where the final u integral evaluates to ln((l — u)/{ 1 + 
u )) + in and contains the imaginary part. 

The real part of the lowest order 3yA energy shift 
AE lo has contributions from the full 0 < p, q < oo quad¬ 
rant. The imaginary part of AElo is a 2-dimensional in¬ 
tegral and comes only from the decay triangle. There was 
no infrared difficulty arising from the loop integrals in¬ 
volving annihilation photons because all six permutations 
of the annihilation photons were considered together. We 
preformed the numerical integrations using the adaptive 
Monte Carlo integration routine Vegas [IT]. Quadruple 
precision was required to control the delicate cancella¬ 
tions among the many parts of the integrands. Our total 
result was 

AE lo = { -0.5126319(20) - 0.30354919(3)«| —. 

77 ( 10 ) 

The real part of AElo is consistent with the earlier eval¬ 
uation of the 3yA contribution to the o-Ps energy at 
0(ma 6 ) HT 1 U 1 lj . and the imaginary part gives the corre¬ 
sponding decay rate through Flo = — 2lm(AELo) with 
a numerical value that agrees with the Ore and Powell re¬ 
sult Flo = t ^( 71 ' 2 — 9 )ma 6 [45]. These checks show that 
our method of calculation is working properly and also 
provide an independent check of the 3 7 A energy shift. 

The one-loop corrections to the lowest order 3 ")A di¬ 
agram are shown in Fig. [2] These terms all give energy 
contributions of 0(ma 7 ). They were evaluated one by 
one as different techniques were required for the vari¬ 
ous contributions. Feynman gauge was used through¬ 
out. The self energy (Fig. [2^,) and outer vertex (Fig. [2jo) 
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FIG. 2: The seven types of one-loop radiative corrections 
in the 3yA channel. Each of these diagrams represents the 
full set of permutations of the annihilation photons and ad¬ 
ditional contributions coming from the various places where 
the correction could act. 


were fairly straightforward. We used known expressions 
in terms of Feynman parameters for the self energy and 
vertex functions [3S1 05; and calculated the imaginary 
parts as integrals over the decay triangle and real parts 
as integrals over the full pq quadrant. The results are 
shown in Table [II] A photon mass A was introduced to 
allow for mass-shell renormalization of the self energy 
and vertex parts and to regulate the binding singularity 
in the ladder graph. All dependence on A vanished in the 
net contribution. 

The inner vertex contribution of Fig. was not 
straightforward. The parametric expression for the ver¬ 
tex function contains a denominator factor that is a com¬ 
plicated quadratic function of po and qo that can’t be 
easily integrated using the residue theorem. Instead, we 
derived an expression for an altered vertex function hav¬ 
ing the same UV behavior (to allow renormalization) but 
with no dependence on qo in the denominator. The dif¬ 
ference is UV-finite, and we evaluated the integral over 
the loop momentum k non-covariantly using our residue 
theorem routine for dk 0 while d 3 k was done as part of 
the overall numerical integration. We found the region 
of small p and q to be particularly challenging for our 
numerical integration routine. We broke up the integra¬ 
tion region into a number of parts near p = q = 0 and 
evaluated them individually. The total integral was the 
result of an extrapolation. 

The double vertex (Fig. [2jl) was straightforward with 
no renormalization nor IR regularization required. We 
performed the loop integral as dko using the residue 
theorem routine and d 3 k numerically. Again an extrap- 
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olation was necessary to achieve an acceptable result for 
small p and q. 

The ladder graph (Fig. [2^) is UV safe but contains a 
binding singularity that required special treatment. The 
IR divergence is regulated by the non-zero photon mass 
we have assumed throughout. We separated off the most 
singular part: 2Ib(oi/'k)A.Elo where the binding integral 


loop (to facilitate implementation of the UV cancella¬ 
tion) seems out of bounds because po and qo would get 
hopelessly mixed up, while a parametric evaluation of the 
graph as a whole would give an expression from which the 
real and imaginary parts would not be easily disentan¬ 
gled. The key to the evaluation of this graph is the gauge 
identity satisfied by the light-by-light tensor )53] 


/ W 4 k — 1 t r 

S?5W = a + 1 " a - 1 + °< a > (11) 

with D(k) = {k 2 — \ 2 ){(k+n) 2 — l)((k—n) 2 — 1). The7r/A 
term is the binding or threshold singularity in this ap¬ 
proach to the problem. Its significance is that the ladder 
correction contains the basic Coulomb binding that holds 
the atom together in the first place. The 7 t/A must be 
discarded, as retaining it would amount to double count¬ 
ing [47] , To make the point in another way, exactly this 
7 t/A would cancel in the matching procedure if the cal¬ 
culation had been done using an effective non relativistic 
field theory formalism 48j [49]. The (In A — 1) multi¬ 
plies an expression proportional to the lowest order 3y A 
contribution (times a/ir). The remainder of the ladder 
contribution is IR finite and was dealt with in the usual 
way. 

The light-by-light graph (Fig. [2^) promised to be a 
challenge. We could not just integrate the loop momen¬ 
tum non-covariantly using our residue theorem routine 
because the light-by-light loop has a UV divergence as 
seen by naive power counting (although this divergence 
vanishes when all permutations are included). Further¬ 
more, the graph as a whole has an apparent UV diver¬ 
gent by the same power counting argument. On the 
other hand, a parametric evaluation of the light-by-light 


TABLE II: Numerical results for the various contributions 
to the energy at 0{ma 7 ) due to one-loop corrections to the 
3yA process. Contributions to I are shown where A E = 
(ma 7 /n 3 )I. The IR singular terms are proportional to Ilo, 
where AElo = {ma 6 /tv 2 )Ilo- The real and “calculated” 
imaginary contributions were obtained here via numerical in¬ 
tegration. The “analytic” imaginary results come from prior 
work in which the o-Ps decay rate contributions were calcu¬ 
lated exactly [501 ;51 or semi-analytically (in terms of two- 
dimensional integrals |52| ). 


Term 

In A Ilo 

Real 

Imaginary 

(calculated) 

Imaginary 

(analytic) 

SE 

4 

-4.48169(23) 

-1.452474(3) 

-1.452478 

OV 

-4 

4.47981(26) 

0.312310(1) 

0.312310 

IV 

-2 

2.28614(44) 

0.558326(2) 

0.558325 

DV 

0 

-1.35555(59) 

1.082948(4) 

1.082944 

LAD 

2 

1.22037(48) 

2.374291(7) 

2.374287 

LbyL 

0 

1.05835(39) 

0.247109(8) 

0.247106 

VP 

0 

-0.585790(2) 

0 

0 

total 

0 

2.62164(103) 

3.122510(12) 

3.122494 


n 


Ml >M2 7 M3 A 


_ _ pK, 


d p\ 


n 




( 12 ) 


where P is the 4-momentum of the incoming photon 
and k is the corresponding Lorentz index. The gauge 
identity is obtained by differentiating the condition of 
gauge invariance as applied to the light-by-light tensor 
P K n M1)At2)/l3i(S = 0. The terms on the right hand side of 


(12) are each well-behaved in the UV and were treated 


by non-covariant integration of the fermion loop momen¬ 
tum. 

While the light-by-light graph contains a virtual an¬ 
nihilation to a single photon, it was not included in the 
earlier evaluation [541 of one-photon-annihilation contri¬ 
butions [54] . 

The vacuum polarization contribution (Fig. [2jr) is 
purely real because the vacuum polarization-corrected 
photon does not have the usual pole at p 2 = 0. We used 
a standard parametrization of the vacuum polarization 
function [3S| and performed the integrals over p and q in 
two ways: by Feynman parameters and non-covariantly 
using our residue theorem routine. The results of the two 
approaches were consistent. 

The results for the various contributions are summa¬ 
rized in Table [TT] The net IR divergence vanishes, as it 
must. The imaginary part of each contribution is in ac¬ 
cord with the corresponding result of prior work. The 
net result of all 0(a ) corrections in the three-photon- 
annihilation channel is 


A E = {2.6216(11) + 3.122510(12) i} 


(13) 


The energy level shift implied by this correction is 


TfloJ Sf r\Sc 1 
A E = 2.6216(11)—— . 

7 t 6 n 6 


(14) 


for states with principal quantum number n, orbital an¬ 
gular momentum £, and total spin s. For o-Ps in the 
ground state this amounts to 11.5 kHz. On comparing 


(10) and (13) we see that these corrections involve siz¬ 


able factors: 


RcAU 3 7j 4 

ImA£ 37 A 


—0.5126 j 1 — 5.114(2) “l 7 "" 6 , (15a) 

1 7T J 7T 2 

0.3035{ 1 - 10.2866-}^-. (15b) 

1 7T J TT Z 


The correction to the imaginary part is easily seen in o-Ps 
lifetime measurements [55] . The 3yA energy correction is 
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relevant to all classes of positronium spectroscopy: 1S-2S, 
hfs, and n = 2 fine structure. It should certainly be seen 
if suggested improvements in the 1S-2S measurements 
can be realized. 
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